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Abstract. Isolated Herbig Ae stars can be divided into two groups (Meeus et al. 1200 J) : those with an almost 
flat spectral energy distribution in the mid-infrared ('group I'), and those with a strong decline towards the far- 
infrared ('group IP). In this paper we show that the group I vs. II distinction can be understood as arising from 
flaring vs. self-shadowed disks. We show that these two types of disks are natural solutions of the 2-D radiation- 
hydrostatic structure equations. Disks with high optical depth turn out to be flaring and have a strong far-IR 
emission, while disks with an optical depth below a certain threshold drop into the shadow of their own puffed-up 
inner rim and are weak in the far-IR. In spite of not having a directly irradiated surface layer, self-shadowed 
disks still display dust features in emission, in agreement with observations of group II sources. We propose an 
evolutionary scenario in which a disk starts out with a flaring shape (group I source), and then goes through the 
process of grain growth, causing the optical depth of the disk to drop and the disk to become self-shadowed (group 
II source). We show that this scenario predicts that the (sub-) millimeter slope of the disk changes from steep 
(small grains) to Rayleigh- Jeans-like (large grains) in the early stages of evolution, so that all group II sources 
are expected to have Rayleigh- Jeans-like slopes, while some group I sources may still have steep (sub-)millimeter 
slopes. 
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1. Introduction 

Herbig Ae/Be stars are thought to be the intermediate 
mass counterparts of T Tauri stars (see e.g. Waters & 
Waelkens Il998h . They have a strong infrared (IR) excess 
arising from warm circumstellar dust. Similar to T Tauri 
stars, this dust is believed to reside in a circumstellar disk. 
Although there is mounting evidence for the disk-like dis- 
tribution of this material, this issue is still not completely 
settled. More importantly, the physics and geometry of 
such disks are not yet fully understood. 

A number of He rbig Ae/Be stars have been studied by 
Meeus et al. (,20(h]). They present ISO spectra combined 
with ground-based photometry of 13 Herbig Ae/Be stars, 
and convincingly show that they can be classified into two 
main groups, one of which can be subdivided into two 
more subgroups. The main division (group I vs. II) distin- 
guishes the sources on the basis of the shape of the over- 
all spectral energy distribution (SED). Group I sources 
have a relatively strong far-IR flux, which is energetically 
comparable with the flux in the near-IR. Group II sources 
show a similar near-IR excess as group I sources, but their 
flux falls off strongly towards the far-IR. The fraction of 



energy reprocessed by circumstellar dust is typically 30- 
50% in group I sources and about 15% to 30% in group 
II sources. Examples of these two different kinds of SEDs 
are shown in Fig. 

In the paper of Meeus et al. it was speculated that the 
distinction between group I sources and group II sources 
might be explained qualitatively by the disk having a flar- 
ing geometry or not. A flaring disk captures and repro- 
cesses more stellar radiation at large radii than flat disks 
do, and therefore naturally has a stronger far-infrared ex- 
cess. However, it remained unclear what physical mech- 
anism could cause some disks to be flaring and others 
to be flat. In fact, mode ls of passive circum stellar disks 
(e.g. D'Alessio et al. ll998l Chia ng fc G oldreich ll997l CG97: 
Dullemond, Dominik & Natta l200ll DDNOll have so far 
shown that these disks consistently have a flaring geome- 
try. 

A possible mechanism that could flatten a disk is dust 
settling. Due to the gravity excerted by the central star, 
dust grains can slowly drift towards the midplane of the 
disk while the gas remains behind. Since the stellar radi- 
ation is absorbed almost exclusively by dust, the photo- 
sphere of the disk flattens as dust settling proce eds. This 
mechanism has been proposed by Chiang et al. l)200l|) as 
an explanation for the mid- and far-IR SEDs of the Herbig 
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Fig. 1. Two examples of SEDs of Her big Ae stars of the 
sample of Meeus et al. showing the qualitative difference 
in the SEDs of group I sources and group II sources. 

stars MWC 480 and HD36121. They found that the SEDs 
of these sources could only be fitted with a CG97 type disk 
if the vertical geometric thickness was artificially reduced 
by a factor of about 3 to 5. They argue that this suggests 
that dust settling is the mechanism behind the group II 
type SEDs. 

In this paper we will demonstrate a more direct, 
though possibly related, way of causing a disk to have 
a group II type SED. It is the effect of self-shadowing. 
If the disk's geometric thickness H/R is a monotonic in- 
creasing function of radius R, then the disk is flaring, and 
can capture stellar radiation at all radii. If, on the other 
hand, the ratio H/R becomes smaller as one goes to larger 
radii, then the disk becomes self-shadowed: the inner disk 
regions cast a shadow over the outer disk regions. As the 
shadow deprives the outer regions from most of their ir- 
radiation flux (only indirect irradiation remains possible, 
but is rather weak), the far-IR emission is strongly sup- 
pressed. This works even without the need of a factor of 
3 to 5 reduction in H , as long as the outer parts of the 
disk are in the shadow of the inner parts. Dust settling 
could be a possible mechanism causing such a shadowing 
(Dullemond & Dominik in prep.), but another possibility 
is the vertical optical depth of the disk, as we will show in 
this paper. 

We will show that an important role in the self- 
shadowing is played by the disk's puffed up inner rim. 
The existence of such a structure was inferred from ob- 
servations of Her big Ae/Be stars in the near-IR. Virtually 
all Herbig Ae/Be stars show a prominent thermal bump 
in their SEDs around 3 /xm. Since a simple flaring disk 
model is not capable of reproducing it, several authors in- 
voke additional non-disk com ponents to explain this fea- 
ture (Miroshnichenko et al. [1997 1999; van den Ancker 
et al. l2000|) . But Natta et al. l|200l|i have drawn at- 
tention to the importance of the inner hole in the disk 
caused by dust evaporation, and the emission from the 



inner rim of the dusty part of the disk. It was shown by 
Dullemond, Dominik & Natta l|200li henceforth DDN01) 
that the model of Chiang & Goldreich, extended by a self- 
consistent (though simplified) description of the inner rim, 
could indeed explain the entire SEDs of Herbig Ae stars, in 
particular those of group I. An important feature of these 
models is that the inner rim of the dusty disk is puffed-up 
as a result of its direct exposure to the light of the central 
star. The rim therefore casts a shadow over part of the 
disk, typically from the location of the rim (~ 0.5 AU) 
out to about ~ 5 AU. Only further out the surface height 
of the disk becomes sufficiently high to cause the upper 
layers of the disk to reach out of the shadow wedge cast 
by the inner rim (see Fig. |H1 left). 

While the DDN01 models work best for group I 
sources, they can also fit group II sources (even without 
dust settling). But these fits require an unlikely termina- 
tion of the disk just outs ide of the shadow cast by the 
inner rim (Dominik et al. 120031) . It seems unlikely that 
all group II sources have such a fine-tuned outer radius. 
For some sources it is also known from other obser vation s 
that the disks are much larger (Mannings & Sargent ll997|) . 
Dominik et al. argue that the inferred smallness of the 
flaring part of the disk may in fact be a hint that these 
disks are fully self-shadowed, i.e. that the inner rim shad- 
ows the entire disk. However, no detailed models of such 
self-shadowed disks were given. 

A consistent discription of self-shadowed (non-flaring) 
disks can only be achieved by employing a 2-D (axisym- 
mctric) or 3-D radiative transfer code to compute the tem- 
perature structure of the disk. By also computing the ver- 
tical density structure of the disk simultaneously, it can be 
investigated if, and under what circumstances non-flaring 
or self-shadowed disks can exist. A firs t attempt in this 
direction was done by Dullemond (|2002L henceforth D02) . 
Using a grey 2D radiative transfer code it was shown that 
disks with a high optical depth flare in the way described 
by DDN01 (inner rim, shadow and outer flaring part). 
But when the optical depth of the disk is below a cer- 
tain threshold, the entire disk drops into the shadow of its 
own inner rim (see Fig. [S] right). The resulting SED turns 
out to be very similar to group II sources: a strong near 
IR bump and a weak far-IR flux with a steep slope. In 
the D02 paper it was therefore suggested that the group 
II type of SED is caused by self-shadowing, instead of a 
mere flattening of a flaring disk. 

In spite of these new findings, a number of important 
questions remain to be answered before firm conclusions 
can be drawn. In particular, the assumption of grey opac- 
ities in D02 is very unrealistic. If the disk is dominated by 
small grains (as indicated by the prominent emission fea- 
tures) the opacity is far from grey. It is possible that real 
opacities will change the picture considerably. Moreover, 
with real opacities we can compare the model to the obser- 
vations in more detail, and we can assign realistic masses 
to our disk models where previously we could only speak of 
the optical depth. Therefore, the present paper addresses 
the following issues: 
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— Does the group I /ll explanation in terms of 
flaring/self-shadowed still hold when real opacities are 
used? If so, for which disk parameters do we get flar- 
ing or self-shadowed disks? Does dust grain growth 
play a role in the division between group I and group 
II sources? 

— Will self-shadowed disks provide the 10 micron fea- 
ture in emission, as seen in most group II sources? In 
the simple picture of a CG97 or DDN01 disk, these 
features are always formed in the superheated layer 
created by dust directly illuminated by the star. In a 
self-shadowed disk the disk surface is not directly irra- 
diated. 

— What do the models predict for the (sub-)mm flux and 
slope for group I and group II disks? The grey models 
were entirely optically thick in the (sub)mm, which 
may not be the case when real opacities are used. 

In this paper we present the first self- consistent 2-D 
models of Herbig Ae star disks with realistic dust opacities, 
and we show that indeed the discrepancy between group 
I and group II sources can be understood naturally in 
terms of flaring and self-shadowed (i.e. non-flaring) disks. 
We also show that through the reduction of the optical 
depth of the disk by the growth of grains, a flaring disk 
can be turned into a self-shadowed disk. This suggests a 
natural evolutionary link between group I and group II 
sources. 

The paper is organized as follows: in section [21 we de- 
scribe the modeling procedure and radiative transfer codes 
used. In section [21 we describe the structure of - and the 
SEDs produced by - two sets of models which run from 
large to low optical depth of the outer disk parts. In sec- 
tion 0] we discuss the results in the framework of shadow- 
ing, rim emission, solid state features and grain growth. 

2. Models 

2.1. Modeling procedure 

The equations for the models of this paper are very simi- 
lar to the ones presented in D02, but while in D02 a grey 
opacity is used for the sake of simplicity, we now use sil i- 
cate grains with a size of a = 0.1/im (Draine & Lee|l984). 
In order not to complicate the models unnecessarily we 
ignore the scattering opacity (we set it to zero). An ex- 
tensive discussion of the effects of scattering in models of 
passiv e reprocessing disks is given in Dullemond & Natta 

To compute the disk structure, we proceed as follows. 
We adopt as our coordinate system the polar coordinates 
R and G (where = means the pole and = it/ 2 is 
the equator), and we assume axial symmetry. The disk 
is presumed to be passive, i.e. non-accreting, so that the 
energy balance is entirely determined by the irradiation 
by the central star. This process of irradiation is modeled 
using a 2-D axisymmetric continuum radiative transfer 
code that solves for the dust temperature as a function of 
R and 6. 



The determination of the density structure requires an 
iteration procedure. An input to the model is the surface 
density distribution T,{R). At the start we set up a rea- 
sonable initial guess for the density distribution p(R, 0) 
that is consistent with this S(i?). We can then apply the 
radiative code to find the gas temperature everywhere. 
Using this temperature distribution, a vertical integration 
of the equation of vertical hydrostatic equilibrium yields 
the new 2-D density structure p(R, 0). After iterating the 
entire procedure of radiative transfer and vertical struc- 
ture a couple of times one finds a solution in which both 
p(R, 0) and T(R, 0) are consistent with each other. Using 
this density and temperature structure a ray-tracer can 
now produce the required SEDs and images. 

There is, however, a major hurdle that needs to be 
overcome before the above schetched procedure can work. 
The disks that we wish to model may have tremendous op- 
tical depths. An optical depth of ry = 10 5 is not uncom- 
mon. Most radiative transfer codes cannot handle such 
optical depths: they will either produce noisy or wrong 
results, or they never converge. As was described in D02, 
to solve this problem we use the code RADICAL, which 
is based on the method of Variable Eddington Tensors 
(VET). This method is the only method presently known 
to be able to treat problems with such enormous opti- 
cal depth without requiring excessive excecution time. For 
1+1D type disk models they have been succesfully ap- 
plied many times already (see e.g. Mal bet et al. 199li 
l200ll Dullemond, van Zadelhoff & Natta 2002). For 2-D 
radiative transfer problems they are slo wly c oming in use 
(e.g. Nakazato, Nakamoto & Umemura 2003: D02). 

The problem with the VET algorithm implemented in 
RADICAL is that it sometimes has slight inaccuracies in 
the resulting temperature profiles. These inaccuracies are 
usually not larger than about 5% to 10%, but since the 
SED is proportional to T 4 , such inaccuracies are amplified 
in the outcoming spectrum. For the structure iteration 
these inaccuracies are less dangerous, as the pressure scale 
height of the disk is proportional to y/T, meaning that the 
error in the density structure of the disk is 2.5% to 5%. 
It seems therefore that RADICAL is appropriate for the 
structure iteration procedure, but not for computing the 
temperature profile for the ray-tracer producing the SED 
and images. 

For the final temperature structure determination (af- 
ter convergence of the density structure) we therefore pre- 
fer to use a different method: an impr oved version of the 
algorithm of Bjorkman & Wood (2001), implemented in a 
code called RADMC. As expected, this method suffers from 
large CPU costs at high optical depth, and in many cases 
produces an strong numerical noise on the temperature 
profile in the most optically thick regions. But this is 
not a major problem since these optically thick regions 
(usually near the equator at small radii) are unobserv- 
ablc. And since RADMC is called only once per model 
(while RADICAL has to be called many times during the 
structure iteration), the large CPU cost of RADMC re- 
mains manageable. RADMC has the advantage of being 
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very accurate and reliable when it comes to the emerg- 
ing spectrum. We have found that for some of the models 
RADICAL slightly overpredicts the 5-8 /jm part of the 
SED, while RADMC produces a sharper transition be- 
tween the mid-infrared-bump and the 10-micron feature, 
in agreement with the SEDs predicted by s emi-analytic 
models (Dullemond, Domini k k, N atta l200l[) . and found 
in nature (e.g. Meeus et al. l200l[) . We verified that the 
density structure produced after the RADMC run is vir- 
tually identical to the one produced with the RADICAL 
code. Also, both codes have been independently tested 
against other codes on a number of test cases, albeit of 
m uch les s extreme optical depths than here (Pascucci et 
al. l2004(l . Both the RADICAL results and the RADMC 
results conserve energy to within an error of at most 5%. 

As final note it should be mentioned that the itera- 
tion procedure, in which the hydrostatic equilibrium code 
and the radiative transfer code are alternately applied, 
may in some cases not reach a perfectly converged solu- 
tion. Minor waves may propagate over the solution from 
one iteration to the next and never damp out completely. 
This problem is related to a kno wn instability operating 
in these disks (Dullemond l200Cj|) and may point to such 
waves propagating the disk in reality as well. For disks 
around Herbig Ae/Be stars, like the models presented in 
this paper, these waves are weak and have virtually no ef- 
fect on the SED or images. For T Tauri stars the problem 
is more serious, and a study of the stability of these disks 
using time-dependent heating and cooling calculations is 
required. For this reason we limit ourselves in this paper 
to Herbig Ae/Be stars only. 

2.2. Setup of the models 

The models of this paper simulate a disk around a star of 
M* = 2.5M , i?» = 2R Q and T* = 10000A'. The surface 
density of the disk as a function of radius E(i?) is defined 
to be: 

E(i2) = X (R/Ro) p (1) 

with Rq taken to be i?o = 200AU. The mass of the disk is 
then: 
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Similar equations hold for the dust mass of the disk Mdust 
in relation to the dust surface density E du st • In general one 
has M dust = 0.01 Mdi s k and E dust = 0.01 E disk , because of 
the gas-to-dust ratio of 100. 

We take the inner radius to be at 0.5AU, which is 
approximately the dust evaporation radius, and we put 
the outer radius at 200AU. We allow the disk to extend a 
bit further than 200AU, but with a very steep powerlaw 
index for the surface density E cx R~ 12 , so that effectively 
the disk has a is 200 AU outer radius, but does not end 
too abruptly there. 

In this paper we present three series of models: 



— Series A: Disks with equal mass (Mdisk = 0.01M©, 
i.e. M dust = O.OOO1M ), but different distribution of 
E(i?) (power law index between p = — 1 ■ • • — 4). The 
grain size is 0.1/im. 

— Series B: Disks with varying mass (between Mdisk = 
0.1---10- 6 M Q , i.e. Md^t = lO- 3 ---lO- 8 M ), but 
with the same distribution of this mass as a function 
of radius (p = —1.5). The grain size is again 0.1/xm. 

— Series BL: Like the B series, but the reduction in 
mass of the disk is now a reduction only in the small 
grain dust mass (0.1/im). As the mass in small grains 
is reduced, the removed dust mass is converted into a 
layer of large (2 mm) grains at the equatorial plane, 
so that the total dust mass (0.1/im and 1 mm grains) 
remains the same (Mj ust = 0.001 M Q ). In this way we 
simulate the process of grain growth, and the subse- 
quent settling of the large grains to the midplane. The 
small grains are assumed not to have settled. 

The parameters of all these models are listed in table ^ 

The A series allows us to see whether self-shadowed 
disks can result from redistributing matter in the disk 
from the outside to the inside. The B series shows whether 
an overall reduction of the optical depth (i.e. mass) can 
cause self-shadowing. The BL series is like the B series, 
but simulates grain grow th and settling (see D'Alessio, 
Calvet & Hartmann l200ll for an earlier study of the effect 
of grain growth on disk models). 

3. Results 

3.1. Structure of the disk 

To show the structure of the computed models, we use con- 
tours of temperature and density in a plot in which the 
2-dimensional spherical coordinates log(-R) and ir/2 — Q 
are on the x- and y-axes respectively. The midplane is at 
the bottom of the figure (6 = 7r/2). These plots show 
the entire disk structure in a clear way. Lines at con- 
stant 7r/2 — (radial rays) are horizontal lines in this 
plot. Flaring and non- flaring disks can therefore easily be 
distinguished by plotting the t = 1 location along these 
lines: flaring disks will be indicated by rising curves while 
non- flaring curves will show flat curves. The resulting tem- 
perature and density structures for models Al and A4 are 
shown in Fig. [2] in the way explained above. In Fig. [3] the 
same models are shown, but this time only the density 
contours. Overplotted are the r = 1 surfaces at 0.55/im 
and 3/im. In these figures one can clearly see the puffed-up 
structure of the inner rim in both models, in agreement 
with the DDN01 model. But at large radii only model 
Al has rising density contours, while model A4 has con- 
stant or even declining density contours. This difference 
is typical of the the difference between flaring and self- 
shadowed disks. The distinction is more clear in the shape 
of the r = 1 surfaces: for the flaring model (Al) they move 
upward toward larger radius, while for the self-shadowed 
model (A4) they remain at constant ir/2 — 8. It is inter- 
esting to note that for the flaring model (Al) the upward 
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Table 1. The parameters of the three series of models 
presented in this paper, p is the power law index for the 
surface density <x R p ). Md us t, small is the mass in 

small (0.1 /im) grains. These grains are evenly distributed 
with the gas in a constant gas-to-dust mass ratio. For the 
A and B models, this ratio is 100. For the BL models 
this ratio is 100,1000,10000 etc., but still constant over 
the disk. The -Md us t,big is the mass in big (2 mm) grains 
which are assumed to be located in a thin midplane layer. 
This parameter is only non-zero for the BL series, and it 
is taken such that Afdust,smaii + Mdust,big = 10~ 3 M Q . In 
this way, the BL series simulates a process of converting 
small grains evenly distributed in the disk into big grains 
located at the midplane. Mdi s k is the total mass of the disk 
(dust+gas) as calculated from the three disk parameters. 
It is therefore not a model parameter. Note that the Bl 
and BL1 models are identical. 



movement of the r = 1 starts around about 3 AU, while 
the curve remains approximately constant below 3 AU. 
This is the shadow of the inner rim that covers the inner- 
most (< 3 AU) part of the disk. This phenomenon was 
predicted by DDN01, and is confirmed here. 

The temperature contours of Fig. [21 are somewhat 
harder to interpret. Clearly the temperature at the disk 
midplane is lower than the temperature of grains above 
the disk's surface. This is the reason the disk has dust 
features in emission. The 'kink' in the temperature profile 
roughly follows the t = 1 surface, since r = 1 defines the 
surface where the direct stellar radiation is absorbed. 

Many of the qualitative structual features of these 
models are similar to the models of D02, and we refer 
to that paper for a deeper discussion and interpretation 
of the temperature and density contours shown here. 

The models from series B shows similar results. Fig. 0] 
shows the stucture plots of models Bl and B5. Model Bl 
has the same kind of flaring shape as model Al. As one 
goes toward lower mass (B5) the disk tends to become 
more self-shadowed, resembling model A4, but with a less 
high inner rim. This shows that decreasing the mass of the 
disk can have a similar effect as changing the power law 



slope of the disk. But the effect is less abrupt. Model B4, 
for example, is intermediate between a flaring and a self- 
shadowed disk, and even model B5 is less self-shadowed 
than model A4 (compare the left panels of Figs. Eland 0J. 
This is because in the B series both the height of the inner 
rim as well as the height of the surface of the disk behind 
it get reduced. The reduction of the disk surface height 
is faster, and that is the reason why from Bl to B6 the 
flaring disk turns into a self-shadowed one. 

We do not show density and/or temperature structure 
figures for the BL series, since the global structure of the 
BL disk models is the same as for the B disk models. The 
only real difference is the presence of the passive midplane 
layer of large grains in the BL models. These grains do not 
affect the temperature structure of the disk. The fact that 
the gas mass of the BL disks are all identical while for 
the B series they decrease from Bl to B6 (see table ^1 
has no effect on the structure, because the gas-to-dust ra- 
tio drops out of the equations for passive disks, as long 
as it remains a global constant of the model. Therefore, 
the temperature structure and, except for a global fac- 
tor, the density structure of the B and BL disk models 
are the same. In the modeling procedure for the BL series 
we in fact use the density structure of the corresponding 
B model, add the midplane layer, do one more radiative 
transfer run with the RADMC code to compute the mid- 
plane grain temperature and finally create the SEDs for 
the BL series. 

3.2. The SEDs of the disks 

In Fig. the spectral energy distributions of the models 
of series A are shown. The SED varies from one that has a 
strong far-IR flux to an SED that is dominated by near-IR 
flux. The models with a strongly flaring outer disk (models 
A1,A2) are the ones that have a strong far-IR flux. The 
radiation of the central star is thereby captured by the 
disk at large radii, and reprocessed into far-IR emission. 
In particular in model Al, most of the mass is in the outer 
regions of the disk. The flaring part has a large covering 
fraction with respect to the central star, and is therefore 
bright. Moreover, the low density in the inner rim region 
weakens the shadow effect, exposes a larger part of the 
outer disk to stellar radiation and flaring, and reduces the 
near-IR flux. 

The models with a self-shadowed geometry (model 
A3,A4) turn out to have a low far-IR flux (Fig. A 
significant fraction of the disk mass resides in the inner 
regions, close to the rim. The inner rim is optically thick 
to large heights and casts a strong shadow over the disk. 
The mass in the outer regions of the disk is so low that 
the disk surface stays below the shadow cast by the inner 
rim. The entire disk is non-flaring (it receives no direct 
radiation from the star) and falls into the shadow, with 
indirect radiation from the upper parts of the inner rim as 
the only remaining heating source. Nevertheless, the sili- 
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Fig. 2. The temperature and density structure as a fuction of spherical radius r and polar angle 8 (tt/2 is the equatorial 
plane, i.e. the bottom of the figure) for model Al (left) and A4 (right). The grey-scale contours (accentuated by the 
dotted contour lines) are density. These are spaced logarithmically in steps of a factor 3.3. The solid contours are 
temperature. The temperature contours are spaced 50 K apart. The small tick mark on one of these contours tags the 
200 K contour. 




1 10 100 1 10 100 

R [AU] R [AU] 

Fig. 3. As Fig.|2 but instead of temperature contours, the figure shows the surface of ro.55^ m , r adiai = 1 (solid), and 
T^nm, radial = 1 (dashed). The shape of these curves clearly show that the left model is flared while the right model is 
self-shadowed. In the flared model it is also seen that the flaring starts for real beyond about 8 AU; between 0.5 and 
8 AU the curve stays at almost constant 6, typical for the shadowed region of the disk. 



cate feature at 10 /xm is still in emission, albeit weaker in 
the A4 model than in the other models. 

The SEDs of the B series are shown in Fig.El The first 
five models are still reasonably optically thick, and are 
representative for Herbig Ae/Be stars. Models B6 starts 
to become optically thin even at near-IR wavelengths, 
and therefore represents a transition toward optically thin 
disks such as the debris disks often found around Vega- 
type stars. It is clear that the flaring disk models (Bl to 
B3) still have a reasonably strong far-IR flux while the 
self-shadowed models (B4 to B6) have a weak far-IR flux. 
This is consistent with the group I/II distinction. Also, all 
models have a strong 10 /i silicate feature in emission, even 
the self-shadowed models (B4 to B6), which is consistent 
with observations. 



Interestingly, the infrared flux at different wavelengths 
decreases sequentially from model Bl to B6. First the mm 
flux diminishes, then the far-IR, followed by the mid-IR 
and finally the near-IR. The near- to far-IR part of this 
trend seems to be consistent with the observed differences 
between group I and group II sources. However, the (sub- 
)mm fluxes of the 'group II models' (B4 to B6) are much 
less than those of real group II sources. In fact, the (sub- 
)mm fluxes of group II sources are in reality similar to 
those of group I sources. Therefore, the group I/II dis- 
tinction can not be explained by merely a distinction in 
disk mass. 

The BL series is similar to the B series but the removed 
mass in small grains is put into large grains located at 
the disk's midplane. In Fig. the SEDs of these models 
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Fig. 4. As Fig. 01 but this time for model Bl (left) and B5 (right). 




Fig. 5. The SEDs of models A1...A4 plotted over each 
other. The SEDs are at inclination i = 45°. 



are shown. The near- to far-IR behavior of the BL series 
is the same as the B series. But in the BL models the 
midplane layer of large grains keeps the (sub-)mm flux 
to a certain level, even when the far-IR flux drops from 
Bl to B6. Therefore, from these results one expects no big 
systematic differences in the (sub-)mm flux levels between 
group I and group II sources, which is consistent with 
observations. 

Another striking feature of the BL series is that model 
BL1 (which only has small grains) has a much steeper 
(sub-)mm slope than the BL2 to BL6 models. If grain 
growth and settling governs the transition from group I 
to group II, we would therefore predict that some (though 
not all) group I objects have steep (sub-)mm slopes, while 
all group II sources have less steep (i.e. more Rayleigh- 
Jeans-like) slopes. This is indeed consistent with observa- 
tions: group I disks sometimes have steep slopes, indicat- 
ing that they are dominated by small grains, while group 
II disks generally have much more shallow slopes indica- 



Fig. 6. The SEDs of models B1...B6 plotted over each 
other. The SEDs are at inclination i = 45°. 



tive of large grains (Bouwman et al. l2000t Acke et al. in 
prep.). 

One thing that is not very well in agreement with ob- 
servations is that in the B and BL series, as one goes from 
flaring to self-shadowed, also the flux around 10 fim is 
somewhat suppressed. The ratio of 10 /zm to 40 /im flux 
does not change too much. In reality there exist many 
group II sources with a much larger flux ratio between 10 
/im and the far-IR than these models predict. 

Nevertheless, the BL series, as a model of the origin of 
the group I versus group II distinction, seems to be reason- 
ably well in agreement with observations, and may point 
toward an evolutionary link between these two groups of 
Herbig Ae/Be stars. 

4. Discussion 

4.1. Flaring versus self-shadowed 

We have found two main types of solutions to the coupled 
2-D continuum radiative transfer and vertical structure 
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Fig. 7. The SEDs of models BL1...BL6 plotted over each 
other. The SEDs are at inclination i = 45°. 



problem for passive irradiated circumstellar disks around 
Herbig Ac/Be stars. Some disks are flared from about 3 
AU outwards, and other disks are entirely self-shadowed. 
A pictographic representation of the two main types of 
solutions is shown in Fig. |SJ The SEDs from these two 
types of models are very reminiscent of the two types of 
SEDs observed from Herbig Ae/Be stars: those with a 
strong far-IR flux and those with a weak far-IR flux. 

It is not uniquely determined which parameter 
switches the disk from one mode into the other. Self- 
shadowing is driven by two main factors: there must be 
enough matter in the inner rim to produce a shadow, and 
there must be too little mass at large radii to keep the 
disk's surface above the shadow. The first condition is al- 
most always met: already for very little matter near the 
inner edge the disk becomes optically thick, and a shadow 
is cast. In fact, if so little matter is present that no shadow 
would be cast, then the disk would presumably be globally 
optically thin, and the system would not have been clas- 
sified as a Herbig Ae/Be star in the first place (it would 
more resemble a Vega- type star) . The second condition is 
not always met: as our results show, there are situations in 
which the disk flares, while in other situations it becomes 
self-shadowed. This depends on how high the surface den- 
sity (more accurately: opacity) is at those large radii. If 
the surface density exceeds a certain threshold, the disk 
starts to flare. This threshold, however, depends on radius 
and on the hight of the shadow. 

If we view these results in terms of the model parame- 
ters Mdisk and p, we can see that there are two ways of ar- 
riving at the different solutions. For a given disk mass, the 
surface density powerlaw will determine the surface den- 
sities in the outer disk. Self-shadowed solutions are then 
favored by steep powerlaws, which make the rim higher 
and reduce the surface densities in the outer regions. If 
on the other hand we keep the surface density powerlaw 
fixed, and vary the mass, then self-shadowed disks will be 
favored by low overall disk masses. 



4.2. Dust settling as the cause of self-shadowing 

Another possible way of turning a flaring disk into a self- 
shadowed disk is dust settling. In this paper we show that 
the self-shadowing geometry can be found even without 
invoking dust settling, but it is clear that dust settling 
will take place and may influence the models. For compact 
grains with a size of about 0.1 fxm the settling time scale 
in the outer regions of a 100 AU disk are of the order of 1 
Myr. This means that in the typical life time of a disk the 
dust settles closer to the midplane. Depending on vertical 
mixing processes this may not proceed all the way to very 
thin pancake-like disk, but the geometric thickness of the 
disk can be significantly reduced. If by that that time the 
shadow of the inner rim is still present, such disks would 
naturally become self-shadowed. 

We leave open the question what is more important 
to produce a self-shadowed disk: disk parameters, or the 
effect of dust settling. In a forthcoming paper we will 
model the process of dust settling in detail (Dullemond 
& Dominik in prep.). 

4.3. Shadowing and the silicate feature 

The fact that in T Tauri stars and Herbig Ae/Be stars 
the dust features are generally seen i n emi ssion has been 
explained by earlier models (Calvet Il99ll CG97 etc) as 
arising from the superheated surface layer created by the 
direct irradiation of the disks's surface by the central star. 
Naively one would therefore expect that the self-shadowed 
disks found in this paper would not have emission features, 
as their surfaces do not receive direct stellar radiation. Yet, 
as can be seen in Figs. 03 El and [7| even the self-shadowed 
disks have their silicate feature strongly in emission. There 
are several reasons for this. First of all, stellar photons 
can reach the shadowed region by indirect means (Fig. . 
For instance, photons scattering off the upper parts of the 
puffed-up inner rim may get diverted into the shadowed 
region, and the thermal near-infrared emission from this 
upper part of the rim may also irradiate the shadowed 
region. In fact in the models presented here only the latter 
is at work, since we have neglected scattering. 

Another way by which the shadowed region can be 
heated occurs when the shadowing is imperfect. The A4 
model has an almost perfect shadow (see Fig.|3J), while the 
shadowing in the B5 model is clearly less complete (see 
Fig. 0J . This means that even though one cannot really 
speak of a flaring disk, there is still a tail of the vertical 
density distribution of the disk that still above the shadow 
and may therfore produce a dust emission feature. 

Finally, when the disk has lost so much of it's small 
grain content that it becomes optically thin at 7 /im and 
13 /zm, the silicate feature naturally appears in emission. 
This in fact plays a role in models B5 and B6. 

The indirect heating of the shadowed regions of the 
disk are crucial, not only for the silicate emission fea- 
ture, but also for the far-IR flux. If the shadowing were 
so perfect that no indirect heating were possible, the self- 
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Fig. 8. Pictographic representation of the two disk geometries found. 




Fig. 9. Pictographic representation of how the shadowed 
regions of the disk can be irradiated in an indirect way 
through thermal emission of (and scattering off) the upper 
part of the puffed-up inner rim. Note that the geometry 
is exaggerated in this pictogram. 

shadowed disks outside of the inner rim would be essen- 
tially cold (meaning that they would have the temper- 
ature of the surrounding molecular cloud). They would 
therefore barely produce far-IR flux at all. This would be 
inconsistent with observations of group II sources which 
clearly do show some far-IR flux, albeit much weaker than 
that of group I sources. The indirect heating of these 
self-shadowed regions is responsible for this weak-but-non- 
negligible far-IR flux, and it is reponsible for the fact that 
our self-shadowed disks of the BL-series have SEDs that 
are so similar to existing group II sources. 

4.4. Extended emission: PAH features and scattering 
images 

An important way to study the structure of circumstel- 
lar disks in general and self-shadowing in particular is 
through resolved images. In the sub-mm wavelength re- 
gions, all models discussed in this paper will be extended 
- even though later models of series A will look more com- 
pact because most of the disk mass has been moved to 
the inner regions. A much better test for self-shadowing 
will be images obtained in near and mid-IR PAH bands, 
and images obtained in scattered light. PAH emission re- 
sults from the excitation of PAHs through individual UV 
photons. Scattered light emerges from the disk when stel- 
lar photons in the optical and near-IR are scattered off 
the disk surface into the line of sight. Both processes are 
not modeled in the current paper, so we cannot quantita- 
tively predict observations. However, since both directly 
track the illumination of the disk surface with stellar pho- 



tons (or, on a much lower level, stellar photons scattered 
towards the disk in the upper regions of the rim), we can 
make qualitative predictions. If a disk is self-shadowed, 
both PAH emission and scattered light should be reduced 
dramatically in the outer disk region, much more so than 
in a disk which is still flaring but merely has a reduced 
surface height. In the latter case, the expected reduction 
factor in the intensity of both PAH and scattered light is 
only a few while in fully self-shadowed disks it should be 
orders of magnitude. We note that the models BL4. . . BL6 
are not perfectly self-shadowed - a tenuous tail of the ver- 
tical density distribution of the disk still reaches out of 
the shadow. However, even a small amount of dust set- 
tling would be enough to move the available scatterers 
fully into the shadow. Possible tests for the effect of shad- 
owing in disks therefore include correlating the SED type 
(group I versus group II) with PAH emission and scattered 
light images. In group II sources, the PAH emission should 
be much weaker or even absent - if present it should be 
compact. Scattered light images of group II sources should 
be very dim or even completely dark, except at the inner 
rim. Full studies of these effects are underway, with first 
results indicating that indeed, group II sources are usually 
compact (van Boekel et al, 2003) or weak in PAH emis- 
sion (Acke et al, in preparation). Sources detected so far 
in scatte red light are consi stently group I (e.g . Grady et 
al. l200ll Grady et al. l2000t Danks et al. 120011) . 

5. Conclusions 

In this paper we systematically investigated the 
structure and the SEDs of passive dusty pro- 
toplanetary disks around Herbig Ae/Be stars. 
The models are available in the form of ASCII 
tables on the internet at the following URL: 
http : //www . mpa-garching . mpg . de/PUBLICATIONS/DATA/ 
radtrans/f lareshadow/ 

The main conclusions of this paper are as follows: 

1. The SEDs of Herbig Ae/Be stars can be quite naturally 
understood in terms of a dusty circumstellar disk that 
is passively reprocessing the radiation of the central 
star. A self-consistent 2-D axisymmetric disk model 
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based on 2-D continuum radiative transfer and verti- 
cal hydrostatic equilibrium seems to be a reasonably 
accurate description of such a disk. 

2. Two kinds of solutions are found: disks with a flaring 
geometry longwards of about 3AU and those who are 
fully self-shadowed by the disk's own puffed-up inner 
rim. The flaring disks have SEDs that seem in agree- 
ment with observed SEDs of group I sour ces (in the 
classification scheme of Meeus et al. 12001^1 . while the 
self-shadowed disks have SEDs similar to those ob- 
served from group II sources. These flared and self- 
shadowed disks are natural solutions of the combined 
equations of radiative transfer and hydrostatics. 

3. The mass distribution within the disk is one impor- 
tant factor in determining whether a disk is flared or 
self-shadowed. The total mass of the disk is impor- 
tant as well, but may conflict with the observed fact 
that group II sources are not systematically less mas- 
sive than group I sources. Keeping the mass fixed, but 
growing a large fraction of the grains (by coagulation 
and settling) into a layer of mm size grains at the mid- 
plane does in fact give the right trend from group I to 
II with a reasonable (sub)mm flux. Moreover, this nat- 
urally leads to a less steep (sub)mm slope for group II 
disks which is indeed observed. Therefore, dust grain 
growth and settling can be the driving cause of an 
evolutionary transition turning flaring disks (which ap- 
pear as group I source) to self-shadowed disks (appear- 
ing as group II sources). 

4. All disks (both the flaring and the self-shadowed ones) 
have a 10 /Ltm silicate feature in emission. This is in 
agreement with the known SEDs of group I and group 
II sources. Only when the disk is seen edge-on, the 
silicate feature can turn into absorption. 

5. The models predict a systematic difference in the 
strength and extendedness of PAH emission features 
between group I and group II sources, with the group II 
sources being weaker than the group I sources. A simi- 
lar effect is predicted for resolved images of the disks in 
scattered light: the group II sources being much weaker 
than the group I sources, perhaps even undetectable. 

6. A 2D treatment of radiative transfer is essential in or- 
der to find and treat self-shadowed disk solutions. 
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